The error reported above for dump_load.cpp also occurs for the Python unit test test_dump_load.py only when using parallel Meep with 2 or more MPI processes. Executing these two tests with serial Meep does actually run the test but then fails due to a fields mismatch assertion. Seems there is a bug when two MPI processes try to write to the same file.

I think there is a bug in test_dump_load.py where the reference run and the load-fields run are offset by a single timestep due to this line:

        # Run until t=1 so that fields are initialized.
        sim.run(mp.at_every(5, get_field_point), until=1)

After removing this line as well as geometry, boundary_layers, and symmetries from the Simulation constructor, the test which compares the fields from the two runs at identical times now passes:

    # Tests dumping/loading of fields & structure.                                                                
    def _load_dump_fields(self, single_parallel_file=True):
        resolution = 50
        cell = mp.Vector3(5, 5)
        sources = mp.Source(src=mp.GaussianSource(1, fwidth=1.0),
                            center=mp.Vector3(),
                            component=mp.Ez)
        sim1 = mp.Simulation(resolution=resolution,
                             cell_size=cell,
                             sources=[sources])

        sample_point = mp.Vector3(0.12, -0.29)
        ref_field_points = {}

        def get_ref_field_point(sim):
            p = sim.get_field_point(mp.Ez, sample_point)
            ref_field_points[sim.meep_time()] = p.real

        # First run until t=25 and save structure/fields                                                          
        sim1.run(mp.at_every(5, get_ref_field_point), until=25)

        dump_dirname = os.path.join(self.temp_dir, 'test_load_dump_fields')
        sim1.dump(dump_dirname,
                  dump_structure=False,
                  dump_fields=True,
                  single_parallel_file=single_parallel_file)

        # Then continue running another 25 until t=50                                                             
        sim1.run(mp.at_every(5, get_ref_field_point), until=25)
        print('ref_field_points = ' + str(ref_field_points))

        # Now create a new simulation and try restoring state.                                                    
        sim = mp.Simulation(resolution=resolution,
                            cell_size=cell,
                            sources=[sources])

        sim.init_sim()
        field_points = {}

        def get_field_point(sim):
            p = sim.get_field_point(mp.Ez, sample_point)
            field_points[sim.meep_time()] = p.real

        # Now load the fields.                                                                                    
        sim.load(dump_dirname,
                 load_structure=False,
                 load_fields=True,
                 single_parallel_file=single_parallel_file)
        sim.run(mp.at_every(5, get_field_point), until=25)
        print('field_points = ' + str(field_points))

        for t, v in field_points.items():
            self.assertAlmostEqual(ref_field_points[t], v)

output

Using MPI version 3.1, 1 processes
Saving temp files to dir: /tmp/meepNNflLc
ssRunning test_load_dump_fields
-----------
Initializing structure...
time for choose_chunkdivision = 1.3597e-05 s
Working in 2D dimensions.
Computational cell is 5 x 5 x 0 with resolution 50
time for set_epsilon = 0.0829253 s
-----------
run 0 finished at t = 25.0 (2500 timesteps)
creating fields output file "/tmp/meepNNflLc/test_load_dump_fields/fields.h5" (1)...
Dumped fields to file: /tmp/meepNNflLc/test_load_dump_fields/fields.h5 (True)
run 1 finished at t = 50.0 (5000 timesteps)
ref_field_points = {50.0: 0.17848575860261917, 35.0: -0.09694628417491913, 20.0: 0.12576791644096375, 5.0: -0.32930123805999756, 40.0: -0.12000752240419388, 25.0: -0.24073712527751923, 10.0: 0.1038854569196701, 45.0: 0.05137018859386444, 30.0: 0.24092617630958557, 15.0: 0.07364945858716965}
-----------
Initializing structure...
time for choose_chunkdivision = 5.71e-06 s
Working in 2D dimensions.
Computational cell is 5 x 5 x 0 with resolution 50
time for set_epsilon = 0.0641066 s
-----------
reading fields from file "/tmp/meepNNflLc/test_load_dump_fields/fields.h5" (1)...
Loaded fields from file: /tmp/meepNNflLc/test_load_dump_fields/fields.h5 (True)
on time step 2500 (time=25), 2.94857 s/step
run 0 finished at t = 50.0 (5000 timesteps)
field_points = {50.0: 0.17848575860261917, 35.0: -0.09694628417491913, 40.0: -0.12000752240419388, 25.0: -0.24073712527751923, 45.0: 0.05137018859386444, 30.0: 0.24092617630958557}
.sss
----------------------------------------------------------------------
Ran 6 tests in 0.671s

OK (skipped=5)

In test_dump_load.py, after replacing the line:

sim.run(mp.at_every(5, get_field_point), until=1)

with

sim.init_sim()

The subtest test_load_dump_fields fails because of theAl (aluminum) material used in one of the two Block objects of the geometry:

        geometry = [mp.Block(material=Al, center=mp.Vector3(), size=one_by_one),
                    mp.Block(material=mp.Medium(epsilon=13), center=mp.Vector3(1), size=one_by_one)]

The reason for this test failure is simple: there is a polarization_state object of the fields class which is not being dumped and loaded along with the other field objects as it needs to be:

As a verification that this is indeed the cause of the failing test, replace material=Al with e.g. material=mp.Medium(index=3.2) in the geometry definition above and test_dump_load.py will then pass.

The slight challenge though for dumping/loading the polarization_state is that it is a typedef:

in which the susceptibility object is a generic class:

Dispersive materials such as Al are lorentzian_susceptibility objects which derived from susceptibility:

There are other derived classes of susceptibility besides lorentzian_susceptibility such as multilevel_susceptibility, gyrotropic_susceptibility, etc but for now we can just support checkpointing lorentzian_susceptiblity objects since they are more common than the other varieties.

Codecov Report

Merging #1738 (b8bcf9d) into master (3da5613) will increase coverage by 0.92%.
The diff coverage is 87.50%.

@@            Coverage Diff             @@
##           master    #1738      +/-   ##
==========================================
+ Coverage   73.46%   74.39%   +0.92%     
==========================================
  Files          13       13              
  Lines        4557     4581      +24     
==========================================
+ Hits         3348     3408      +60     
+ Misses       1209     1173      -36     

The fact that you have to use init_sim, and that restoring the checkpoint does not work if you run for a nonzero amount of time first, seems to indicate that there is some state we are forgetting about. It would be good to chase this down.

I think I know why init_sim() is necessary. To understand this involves going through the sequence of allocating the fields arrays that Meep is using. The important point is that the field allocation (i.e., for Ez which is the cause of the segmentation fault reported above) during init_sim() actually happens when calling add_sources:

add_sources in turn calls the member function require_source_components() of the fields class:

require_source_components then calls the lower-level function _require_component for each of the fields required by the particular source component. In the unit test _load_dump_fields, the source component is Ez (which is actually an electric current J_z) and therefore the components in the 2d cell are Ez, Bx, and By. For μ=1 (non-magnetic materials), Meep uses Hx=Bx and Hy=By to avoid unnecessary duplication. Note that in 2d there are two orthogonal sets of polarizations (Ex,Ey,Hz vs. Hx,Hy,Ez) and typically only one set is used in the simulation.

Finally, _require_component calls alloc_f to allocate the fields:

What needs to happen then when dumping and loading the fields is that only those components generated by the particular source need to be saved and allocated. This means that we should save the source component as part of the checkpoint and use it to allocate the correct set of fields when loading the fields similar to how it is already set up.

It seems there is a memory leak in both structure::load and fields:load (introduced in #1715 and this PR, respectively) which is triggered by the C++ unit test dump_test.cpp as reported by the logs of the address sanitizer.